Overall Objective

Load Libraries

library(tidyverse)
library(cowplot)
library(broom)
library(plotly)

Import data

Convert data from ‘wide’ to ‘long’ format

# data
data1 <- data %>%
  gather(Sample,Count,2:49)

# Separate samples by identifiers 
data2 <- data1 %>% 
  separate(Sample, into=c("Sample_ID","Dilution_factor",
                          "Injection","Tech_rep", sep = "_")) %>% 
  select(-`_`)

# Standards
standards1 <- standards %>% 
  gather(Sample,Count,2:25)


standards2 <- standards1 %>% 
  separate(Sample, into=c("Sample_ID","When","Dilution_factor",
                          "Nano_day","Injection","Tech_Rep", sep = "_")) %>% 
  select(-`_`)

Factor the data into categorical variables

# Refactoring Columns for samples
data2$Sample_ID <- as.factor(data2$Sample_ID)
data2$Dilution_factor <- as.numeric(data2$Dilution_factor)
data2$Injection<- as.factor(data2$Injection)
data2$Tech_rep <- as.numeric(data2$Tech_rep)

data2
## # A tibble: 48,000 × 6
##    particle_size Sample_ID Dilution_factor Injection Tech_rep Count
## *          <dbl>    <fctr>           <dbl>    <fctr>    <dbl> <int>
## 1            0.5         1              25         1        0     0
## 2            1.5         1              25         1        0     0
## 3            2.5         1              25         1        0     0
## 4            3.5         1              25         1        0     0
## 5            4.5         1              25         1        0     0
## 6            5.5         1              25         1        0     0
## 7            6.5         1              25         1        0     0
## 8            7.5         1              25         1        0     0
## 9            8.5         1              25         1        0     0
## 10           9.5         1              25         1        0     0
## # ... with 47,990 more rows
# Refactoring COlumns for key
key$Sample_ID <- as.factor(key$Sample_ID)
key$Time <- as.factor(key$Time)
key$Treatment <- as.factor(key$Treatment)
key$Volume <- as.numeric(key$Volume)

key$Treatment <- factor(key$Treatment,levels = c('DMSO','EGF','BPS','BPS_EGF'))

key
## # A tibble: 8 × 4
##   Sample_ID   Time Treatment Volume
##      <fctr> <fctr>    <fctr>  <dbl>
## 1         1     48       BPS    145
## 2         2     48      DMSO    112
## 3         3     48   BPS_EGF    110
## 4         4     48       EGF    120
## 5         5     96      DMSO    105
## 6         6     96   BPS_EGF    112
## 7         7     96       EGF    105
## 8         8     96       BPS    110
# Refactoring columns for standards
standards2$Sample_ID <- as.factor(standards2$Sample_ID)
standards2$When <- as.factor(standards2$When)
standards2$Dilution_factor <- as.numeric(standards2$Dilution_factor)
standards2$Injection <- as.factor(standards2$Injection)
standards2$Nano_day <- as.numeric(standards2$Nano_day)

standards2
## # A tibble: 24,000 × 8
##    particle_size Sample_ID   When Dilution_factor Nano_day Injection
## *          <dbl>    <fctr> <fctr>           <dbl>    <dbl>    <fctr>
## 1            0.5       std  after             125        1         1
## 2            1.5       std  after             125        1         1
## 3            2.5       std  after             125        1         1
## 4            3.5       std  after             125        1         1
## 5            4.5       std  after             125        1         1
## 6            5.5       std  after             125        1         1
## 7            6.5       std  after             125        1         1
## 8            7.5       std  after             125        1         1
## 9            8.5       std  after             125        1         1
## 10           9.5       std  after             125        1         1
## # ... with 23,990 more rows, and 2 more variables: Tech_Rep <chr>,
## #   Count <int>

Back calculate standards

standards2 <- standards2 %>% 
  mutate(True_Count=Dilution_factor*Count)

# Set the correct order of 'categorical factors'
standards2$Nano_day <-  factor(standards2$Nano_day, levels=c('1','2'))
standards2$When <- factor(standards2$When, levels=c('before','after'))

standards2
## # A tibble: 24,000 × 9
##    particle_size Sample_ID   When Dilution_factor Nano_day Injection
##            <dbl>    <fctr> <fctr>           <dbl>   <fctr>    <fctr>
## 1            0.5       std  after             125        1         1
## 2            1.5       std  after             125        1         1
## 3            2.5       std  after             125        1         1
## 4            3.5       std  after             125        1         1
## 5            4.5       std  after             125        1         1
## 6            5.5       std  after             125        1         1
## 7            6.5       std  after             125        1         1
## 8            7.5       std  after             125        1         1
## 9            8.5       std  after             125        1         1
## 10           9.5       std  after             125        1         1
## # ... with 23,990 more rows, and 3 more variables: Tech_Rep <chr>,
## #   Count <int>, True_Count <dbl>

Summarize three technical standard replicates

standards3 <- standards2 %>% 
  group_by(particle_size,Sample_ID,When,Dilution_factor,Nano_day,Injection) %>% 
  summarise( tech_N = length(True_Count),
             tech_mean = mean(True_Count),
             tech_sd = sd(True_Count),
             tech_se = tech_sd/sqrt(tech_N))
standards3
## Source: local data frame [8,000 x 10]
## Groups: particle_size, Sample_ID, When, Dilution_factor, Nano_day [?]
## 
##    particle_size Sample_ID   When Dilution_factor Nano_day Injection
##            <dbl>    <fctr> <fctr>           <dbl>   <fctr>    <fctr>
## 1            0.5       std before             125        1         1
## 2            0.5       std before             125        1         2
## 3            0.5       std before             125        2         1
## 4            0.5       std before             125        2         2
## 5            0.5       std  after             125        1         1
## 6            0.5       std  after             125        1         2
## 7            0.5       std  after             125        2         1
## 8            0.5       std  after             125        2         2
## 9            1.5       std before             125        1         1
## 10           1.5       std before             125        1         2
## # ... with 7,990 more rows, and 4 more variables: tech_N <int>,
## #   tech_mean <dbl>, tech_sd <dbl>, tech_se <dbl>

Summarize standards by injection

standards4 <- standards3 %>% 
  group_by(Nano_day,When,particle_size) %>% 
  summarise( inj_N = length(tech_mean),
             inj_mean = mean(tech_mean),
             inj_sd = sd(tech_mean),
             inj_se = inj_sd/sqrt(inj_N))
standards4
## Source: local data frame [4,000 x 7]
## Groups: Nano_day, When [?]
## 
##    Nano_day   When particle_size inj_N inj_mean inj_sd inj_se
##      <fctr> <fctr>         <dbl> <int>    <dbl>  <dbl>  <dbl>
## 1         1 before           0.5     2        0      0      0
## 2         1 before           1.5     2        0      0      0
## 3         1 before           2.5     2        0      0      0
## 4         1 before           3.5     2        0      0      0
## 5         1 before           4.5     2        0      0      0
## 6         1 before           5.5     2        0      0      0
## 7         1 before           6.5     2        0      0      0
## 8         1 before           7.5     2        0      0      0
## 9         1 before           8.5     2        0      0      0
## 10        1 before           9.5     2        0      0      0
## # ... with 3,990 more rows

Plot before and after plots, facet by experimental day

std_plot <- standards4 %>% 
  ggplot(aes(x = particle_size, y = inj_mean, color=When))+
  geom_line(size=2) + xlim(0,300)+ #line size, x-axis scale
  geom_ribbon(aes(ymin=inj_mean-inj_se, ymax=inj_mean+inj_se),
              alpha=0.4,fill = alpha('grey12', 0.2)) + #error bars
  scale_y_continuous(expand=c(0,0))+ #set bottom of graph
  xlab("Particle Size") + # X axis label
  ylab("\nMean Particle Concentration/ml\n") + # Y axis label
  ggtitle("Nanosight Histogram of\n100nm Standards")+ #title
  labs(color="Condition")+ #Label table title
  facet_wrap(~ Nano_day)

std_plot
## Warning: Removed 1400 rows containing missing values (geom_path).

# ggsave("Standards_histogram_plot.png",
#        height = 5, width = 7, dpi = 300, units= "in")

Standards particle concentrations from each experimental day

standards_df <- standards4 %>% 
  group_by(Nano_day,When) %>% 
  summarise(total=sum(inj_mean))

standards_df
## Source: local data frame [4 x 3]
## Groups: Nano_day [?]
## 
##   Nano_day   When       total
##     <fctr> <fctr>       <dbl>
## 1        1 before 45989600958
## 2        1  after 41407513229
## 3        2 before 51412158250
## 4        2  after 50587138812

Bar graph of standards particle concentrations

standards_bar <- standards_df %>% 
  ggplot(aes(x=Nano_day,y=total,fill=When))+
  geom_col(position="dodge")+
  scale_y_continuous(expand=c(0,0))+ #set bottom of graph
  xlab("Experimental Day") + # X axis label
  ylab("\nMean Particle Concentration/ml\n") + # Y axis label
  ggtitle("Nanosight Histogram of\n100nm Standards")+ #title
  labs(color="When") #Label table title

standards_bar

# ggsave("Standards_bar_plot.png",
#        height = 5, width = 7, dpi = 300, units= "in")

Intraassay variability

Intra.assay_cv <- standards_df %>% 
  group_by(Nano_day) %>% 
  summarise(Intra_Day_N = length(total),
            Intra_Day_mean = mean(total),
            Intra_Day_sd = sd(total),
            Intra_Day_se = Intra_Day_sd/sqrt(Intra_Day_N),
            Intra_Day_cv = Intra_Day_sd/Intra_Day_mean )
Intra.assay_cv
## # A tibble: 2 × 6
##   Nano_day Intra_Day_N Intra_Day_mean Intra_Day_sd Intra_Day_se
##     <fctr>       <int>          <dbl>        <dbl>        <dbl>
## 1        1           2    43698557094   3240025305   2291043865
## 2        2           2    50999648531    583376839    412509719
## # ... with 1 more variables: Intra_Day_cv <dbl>
# # Save as .csv
# write_csv(Intra.assay_cv,"Intra.assay_cv.csv")

Interassay variability

Inter.assay_cv <- Intra.assay_cv %>% 
  summarise(Inter_Day_N = length(Intra_Day_mean),
            Inter_Day_mean = mean(Intra_Day_mean),
            Inter_Day_sd = sd(Intra_Day_mean),
            Inter_Day_se = Inter_Day_sd/sqrt(Inter_Day_N),
            Inter_Day_cv = Inter_Day_sd/Inter_Day_mean )
Inter.assay_cv
## # A tibble: 1 × 5
##   Inter_Day_N Inter_Day_mean Inter_Day_sd Inter_Day_se Inter_Day_cv
##         <int>          <dbl>        <dbl>        <dbl>        <dbl>
## 1           2    47349102812   5162651266   3650545719    0.1090338
# # Save as .csv
# write_csv(Inter.assay_cv,"Inter.assay_cv.csv")

Sample analysis

Back calculate the original concentration of the sample

data2 <- data2 %>% 
  mutate(True_Count = Dilution_factor*Count)
data2
## # A tibble: 48,000 × 7
##    particle_size Sample_ID Dilution_factor Injection Tech_rep Count
##            <dbl>    <fctr>           <dbl>    <fctr>    <dbl> <int>
## 1            0.5         1              25         1        0     0
## 2            1.5         1              25         1        0     0
## 3            2.5         1              25         1        0     0
## 4            3.5         1              25         1        0     0
## 5            4.5         1              25         1        0     0
## 6            5.5         1              25         1        0     0
## 7            6.5         1              25         1        0     0
## 8            7.5         1              25         1        0     0
## 9            8.5         1              25         1        0     0
## 10           9.5         1              25         1        0     0
## # ... with 47,990 more rows, and 1 more variables: True_Count <dbl>

Average three technical readings

data3 <- data2 %>% 
  group_by(particle_size,Sample_ID,Dilution_factor,Injection) %>% 
  summarise( tech_N = length(True_Count),
             tech_mean = mean(True_Count),
             tech_sd = sd(True_Count),
             tech_se = tech_sd/sqrt(tech_N))
data3
## Source: local data frame [16,000 x 8]
## Groups: particle_size, Sample_ID, Dilution_factor [?]
## 
##    particle_size Sample_ID Dilution_factor Injection tech_N tech_mean
##            <dbl>    <fctr>           <dbl>    <fctr>  <int>     <dbl>
## 1            0.5         1              25         1      3         0
## 2            0.5         1              25         2      3         0
## 3            0.5         2              25         1      3         0
## 4            0.5         2              25         2      3         0
## 5            0.5         3              25         1      3         0
## 6            0.5         3              25         2      3         0
## 7            0.5         4              25         1      3         0
## 8            0.5         4              25         2      3         0
## 9            0.5         5              25         1      3         0
## 10           0.5         5              25         2      3         0
## # ... with 15,990 more rows, and 2 more variables: tech_sd <dbl>,
## #   tech_se <dbl>

Summarize samples by injection (average both injections)

data4 <- data3 %>% 
  group_by(particle_size,Sample_ID,Dilution_factor) %>% 
  summarise( inj_N = length(tech_mean),
             inj_mean = mean(tech_mean),
             inj_sd = sd(tech_mean),
             inj_se = inj_sd/sqrt(inj_N))
data4
## Source: local data frame [8,000 x 7]
## Groups: particle_size, Sample_ID [?]
## 
##    particle_size Sample_ID Dilution_factor inj_N inj_mean inj_sd inj_se
##            <dbl>    <fctr>           <dbl> <int>    <dbl>  <dbl>  <dbl>
## 1            0.5         1              25     2        0      0      0
## 2            0.5         2              25     2        0      0      0
## 3            0.5         3              25     2        0      0      0
## 4            0.5         4              25     2        0      0      0
## 5            0.5         5              25     2        0      0      0
## 6            0.5         6              25     2        0      0      0
## 7            0.5         7              25     2        0      0      0
## 8            0.5         8              25     2        0      0      0
## 9            1.5         1              25     2        0      0      0
## 10           1.5         2              25     2        0      0      0
## # ... with 7,990 more rows
# Average technical replicates and merge with key
merge <- left_join(key,data3, by= "Sample_ID")

merge
## # A tibble: 16,000 × 11
##    Sample_ID   Time Treatment Volume particle_size Dilution_factor
##       <fctr> <fctr>    <fctr>  <dbl>         <dbl>           <dbl>
## 1          1     48       BPS    145           0.5              25
## 2          1     48       BPS    145           0.5              25
## 3          1     48       BPS    145           1.5              25
## 4          1     48       BPS    145           1.5              25
## 5          1     48       BPS    145           2.5              25
## 6          1     48       BPS    145           2.5              25
## 7          1     48       BPS    145           3.5              25
## 8          1     48       BPS    145           3.5              25
## 9          1     48       BPS    145           4.5              25
## 10         1     48       BPS    145           4.5              25
## # ... with 15,990 more rows, and 5 more variables: Injection <fctr>,
## #   tech_N <int>, tech_mean <dbl>, tech_sd <dbl>, tech_se <dbl>
# Save as .csv
# write_csv(merge,"Technical_replicate_average.csv")
 
# Average injection replicates and merge with key
merge1 <- left_join(key,data4, by= "Sample_ID")

merge1
## # A tibble: 8,000 × 10
##    Sample_ID   Time Treatment Volume particle_size Dilution_factor inj_N
##       <fctr> <fctr>    <fctr>  <dbl>         <dbl>           <dbl> <int>
## 1          1     48       BPS    145           0.5              25     2
## 2          1     48       BPS    145           1.5              25     2
## 3          1     48       BPS    145           2.5              25     2
## 4          1     48       BPS    145           3.5              25     2
## 5          1     48       BPS    145           4.5              25     2
## 6          1     48       BPS    145           5.5              25     2
## 7          1     48       BPS    145           6.5              25     2
## 8          1     48       BPS    145           7.5              25     2
## 9          1     48       BPS    145           8.5              25     2
## 10         1     48       BPS    145           9.5              25     2
## # ... with 7,990 more rows, and 3 more variables: inj_mean <dbl>,
## #   inj_sd <dbl>, inj_se <dbl>
# #Save as .csv
# write_csv(merge1,"Injection_replicate_average.csv")

Quick visualizations

Graphing all samples

sample_plot <- merge %>%
  ggplot(aes(x=particle_size, y=tech_mean,color=Injection ))+ #plot
  geom_ribbon(aes(ymin=tech_mean-tech_se,
                  ymax=tech_mean+tech_se),
                  alpha=0.2,fill = alpha('grey12', 0.2)) + #error bars
  geom_line(size=2.0, alpha = 0.8) + xlim(0,500)+ #line size, x-axis scale
  scale_y_continuous(expand=c(0,0))+ #set bottom of graph
  xlab("Particle Size") + # X axis label
  ylab("\nMean Particle Concentration/ml\n") + # Y axis label
  ggtitle("Nanosight Histogram of\nhCTBS treated with BPS")+ #title
  labs(color="Injection")+ #Label table title
  facet_grid(Time ~ Treatment)
  # geom_vline(xintercept = 200)+
  # annotate("text", x= 350, y = 1E8, label= "200nm")

sample_plot

# ggsave("Nanosight_Sample_Histogram.png", plot = sample_plot,
#        height = 10, width = 14, dpi = 200, units= "in")

Interactive Plot

ggplotly(sample_plot)

Particle concentration values for each of the samples

merge2 <- merge1 %>% 
  group_by(Time, Treatment, Volume) %>% 
  summarise(particle_conc=sum(inj_mean))
merge2
## Source: local data frame [8 x 4]
## Groups: Time, Treatment [?]
## 
##     Time Treatment Volume particle_conc
##   <fctr>    <fctr>  <dbl>         <dbl>
## 1     48      DMSO    112    5332908929
## 2     48       EGF    120   13941558392
## 3     48       BPS    145    5733223262
## 4     48   BPS_EGF    110    5658288538
## 5     96      DMSO    105    5570460917
## 6     96       EGF    105    4416285746
## 7     96       BPS    110    5557785942
## 8     96   BPS_EGF    112    4699700962

Correct for resuspension volume

merge3 <- merge2 %>% 
  mutate(particle_count = (Volume/1000)*particle_conc, # Create new column with number of particles
         corrected_particle_conc = (particle_conc/1E9)) # Create new column with correct particle concentration
merge3
## Source: local data frame [8 x 6]
## Groups: Time, Treatment [8]
## 
##     Time Treatment Volume particle_conc particle_count
##   <fctr>    <fctr>  <dbl>         <dbl>          <dbl>
## 1     48      DMSO    112    5332908929      597285800
## 2     48       EGF    120   13941558392     1672987007
## 3     48       BPS    145    5733223262      831317373
## 4     48   BPS_EGF    110    5658288538      622411739
## 5     96      DMSO    105    5570460917      584898396
## 6     96       EGF    105    4416285746      463710003
## 7     96       BPS    110    5557785942      611356454
## 8     96   BPS_EGF    112    4699700962      526366508
## # ... with 1 more variables: corrected_particle_conc <dbl>
# Save as .csv
# write_csv(merge3,"Adjusted_particle_concentration.csv")

Barplot

plot1 <- merge3 %>% 
  ggplot(aes(x = Treatment, y = corrected_particle_conc, fill = Treatment)) +
  geom_bar(aes(text = paste("Particle Concentration:",
                            corrected_particle_conc)),
           stat="identity", position = "dodge")+
  xlab("\nTreatment\n") + # X axis label
  ylab("\nMean Vessicle Concentration\n(10^9 particles/ ml)\n") + # Y axis label
  ggtitle("Effect of BPS on Extracellular Vessicle\nRelease of hCTBs\n")+
  scale_y_continuous(breaks = seq(0,14,2),
                     limits = c(0,14),
                     expand = c(0,0))+ # set bottom of graph
  labs(color="Condition")+ # Label table title
  facet_wrap(~Time)
## Warning: Ignoring unknown aesthetics: text
plot1

# ggsave("BPS_treated_hCTBs_48_96_facet_plot.png",
#       height = 8, width = 10, dpi = 600, units= "in")

Interactive Plot

ggplotly(plot1)

Barplot

plot2 <- merge3 %>%
  group_by(Time) %>% 
  ggplot(aes(x = Treatment, y = corrected_particle_conc, fill = Time )) +
  geom_bar(position = "dodge", stat = "identity")+
  xlab("\nTreatment\n") + # X axis label
  ylab("\nMean Vessicle Concentration\n(10^9 particles/ ml)\n") + # Y axis label
  ggtitle("Effect of BPS on Extracellular Vessicle\nRelease of hCTBs\n")+
  scale_y_continuous(breaks = seq(0,14,2),
                     limits = c(0,14),
                     expand = c(0,0))+ # set bottom of graph
  labs(fill= "Time (hr)")

plot2

 # ggsave("BPS_treated_hCTBs_48_96_plot.png",
 #       height = 5, width = 7, dpi = 600, units= "in")

Interactive Plot

ggplotly(plot2)

Statistics

fit <- aov(corrected_particle_conc ~ Time * Treatment ,data=merge3)

tidy(fit)
##             term df    sumsq    meansq
## 1           Time  1 13.57660 13.576598
## 2      Treatment  3 21.35316  7.117719
## 3 Time:Treatment  3 32.29186 10.763954